clear
%% Problem parameters
%% Problem
%   'EqualMinima','Himmelblau','Sixhump','ModifiedRastrigin','Vincent',
%   'Shubert','Composition','IncreasingMinima','Rastrigin','Schaffer'.
ProblemSet.FuncType = 'Shubert';
ProblemSet.dimension = 3;
% ProblemSet.k = [3,3,3];
% ProblemSet.BaseFunc = 3;

[d, x_bound, ObjFunc, OptSet, precision_init, precision, ~, MaximalBudget, FileName]...
    = problem_setting(ProblemSet);

MaximalBudget = 1000000;
iepsilon = 1;
rep=1;
epsilon = [0.1,0.01,0.001,0.0001];
%% Problem
NewRegionNum = 2;
MaxPartitionDepth_d = ceil(log((x_bound(2,:)-x_bound(1,:))./precision_init)...
    ./log(NewRegionNum));
MaxPartitionDepth = sum(MaxPartitionDepth_d)
ClearRadius0 = (x_bound(2,:)-x_bound(1,:))./(NewRegionNum.^(MaxPartitionDepth_d));
% ClearRadius = 2*sqrt(sum(ClearRadius0.^2)); % the maximum distance within a non-partitionable region
ClearRadius = 2*min(ClearRadius0);

Problem.Dimension = d;
Problem.Domain = reshape(x_bound,[1,d*2]);
Problem.Partition = @(Region,SampleSet)Partition_BoxConstraint(Region, NewRegionNum, MaxPartitionDepth, SampleSet);
Problem.Sampling = @(n, region)Sampling_BoxConstraint(n, region, ObjFunc);

%% optimization
% Algorithm parameters
AlgorithmM.n0 = 4;
AlgorithmM.QuantileLevel = 0.3;
AlgorithmM.NewBudget = 200%3;
AlgorithmM.MaxSampleSize = 10; %5*(d+1);
AlgorithmM.radius = ClearRadius;
AlgorithmM.local_search = @(starts,step)LocalSearch_BoxConstraint(ObjFunc,x_bound,starts,step,precision(iepsilon)./2,MaximalBudget);
AlgorithmM.StopCriteria = [1,MaximalBudget];
% AlgorithmM.StopCriteria = [2,size(OptSet.sol,1)];
% AlgorithmM.StopCriteria = [3,5000];

%% STAGE ONE PRS-MMO
time_start = tic;
for i=1:rep
[optima, RegionSet, SampleSet, RegionSampleId, TotalBudget_Stage2] = prsmmo_ls_v3_smallcomparison( Problem, AlgorithmM );
end
toc(time_start)
TotalBudget = size(SampleSet,1)
TotalBudget_Stage2

%% save data
% save(FileName)

%% FIGURES
if rep == 1
    if d==1
        x=x_bound(1,1):(x_bound(2,1)-x_bound(1,1))/1000:x_bound(2,1);
        z=ObjFunc(x');
        z_max=max(z);
        z_min=min(z);
        figure()
        hold on
        plot(x,z)
        for i=1:length(RegionSet)
            plot([RegionSet(i,3), RegionSet(i,3)], [z_min, z_max],'b')
            plot([RegionSet(i,4), RegionSet(i,4)], [z_min, z_max],'b')
        end
        plot(SampleSet(:,2),SampleSet(:,1),'.');
        hold off
    end
    if d == 2
        %% Stage one figures
        figure()
        hold on
        rectangle('Position',[Problem.Domain([1,3]),Problem.Domain([2,4])-Problem.Domain([1,3])])
        scatter(SampleSet(:,2),SampleSet(:,3),2,SampleSet(:,1),'filled');
        for i=1:length(RegionSet)
            rectangle('Position',[RegionSet(i,[3,5]),RegionSet(i,[4,6])-RegionSet(i,[3,5])])
        end
        plot(optima(:,3),optima(:,4),'x','MarkerSize',6,'LineWidth',1.5,'MarkerEdgeColor' ,[0.86,0.34,0.07]);
        hold off
        colorbar
        title('','Interpreter','latex','String',['$N=',num2str(size(SampleSet,1)),'$'])
        xlabel('','Interpreter','latex','String','$x_1$')
        ylabel('','Interpreter','latex','String','$x_2$')
        set(gca,'Fontname','Times New Roman','FontSize',16);
    end
    if d==3
        %% Stage one figures
        figure()
        scatter3(optima(:,3),optima(:,4),optima(:,5),4,optima(:,2),'filled');
        xlabel('','Interpreter','latex','String','$x_1$')
        ylabel('','Interpreter','latex','String','$x_2$')
        title('local search')
        set(gca,'Fontname','Times New Roman','FontSize',16);
%         figure()
%         scatter3(OptSet.sol(:,1),OptSet.sol(:,2),OptSet.sol(:,3),4,'filled');
%         xlabel('','Interpreter','latex','String','$x_1$')
%         ylabel('','Interpreter','latex','String','$x_2$')
%         title('real optima')
%         set(gca,'Fontname','Times New Roman','FontSize',16);
    end
end